Analysis of time-series data using singularities

ABSTRACT

A method for acquiring and analyzing time-series data using singularities is described. This method allows for the analysis of data over a wide spectrum of frequencies. Once the data is acquired in an oil field, singularities of the data are extracted; and the extracted singularities are utilized to interpret the formation properties related to the data.

BACKGROUND

This disclosure relates to a method for analysis and interpretation oftime-series data. More specifically, this disclosure relates to a methodfor analysis and interpretation of time-series data by determining thesingularities of the data. One particular use is in analysis of data,such as acoustic waveforms, from a borehole in a logging operation.

Oil can be found in certain geologic formations at varying depths in theearth's crust. The oil is usually found trapped in layers of porousformations which lies beneath a dome shaped or folded layer of somenon-porous rock such as shale. In other instances, oil is trapped at afault, or break in the layers of the crust.

In the dome and folded formations, natural gas is usually situatedbetween the non-porous layer and the oil described above. The layerbeneath the oil is porous and usually is saturated with water. As iswell known, oil may be collected by drilling a well. This drillingpunctures the reservoir layer and if the penetration is shallow, onlynatural gas will be collected, if the penetration goes too far, then toomuch water will be produced. Therefore, it is highly desirable todetermine the formation properties at a given depth in the borehole.

Also, the composition of the formation surrounding a borehole may be ofinterest. Depending on the formation, it may be evident that anotherarea nearby will provide a better return when drilled. Data acquiredduring logging operations are used to give insight to these materials.

Accurate analysis of acoustic data gathered from a borehole in an oilwell is challenging and complex. The field of sonic logging of boreholesinvolves making acoustic measurements in the borehole at a wide range offrequencies. Acoustic data is generally collected using at least onetransmitter and one receiver. There are several different configurationsin which borehole acoustic data can be collected. These includecross-well imaging (transmitter in one borehole and receiver inanother), borehole seismic (transmitter on the surface and receiver in aborehole), and single-well imaging (transmitter and receiver in the sameborehole). Gathering, separating, and analyzing this acoustic data hassignificant practical applications including fracture identification,compartmentalization, and composition determination.

In order to collect sonic logging data, an acoustic source radiates acompressional wave pulse, which propagates into the surroundingformation. As this pulse enters the formation, it excites compressionaland shear waves in the formation. These waves produce head waves in theborehole fluid that may be measured. Also, in their propagation throughthe formation, the compressional and shear waves encounter interferencethat results in returning energy back towards the borehole where thereceiver may be located. The acoustic responses include head waves andguided borehole waves and Stoneley waves. All of these waves arerecorded by the receiver.

The data gathered by the receiver can provide a large amount ofinformation that is highly valuable for the exploration and developmentof hydrocarbon resources. However, this data must first be accuratelyanalyzed and interpreted. The data retrieved in sonic well logging areextremely complicated because various wave components overlap in timeand in frequency or in both domains simultaneously. Unfortunately,standard Fourier transform techniques are not able to discriminatecomponents that overlap in frequency domain, making it difficult toextract information in this case. In order to work around this weakness,conventional systems typically separate frequency spectra to provide lowfrequency and high frequency analysis. This analysis is useful, butvaluable information can be lost due to the overlapping time andfrequency components.

Accordingly, it would be desirable to envision a method that would allowanalysis of the time series data across a wide range of frequencies,also referred to as a multi-scale analysis. This would allow for a morethorough and informative result that in the instant case could help inanalyzing the formation properties surrounding the borehole.

The difficulties and limitations suggested in the preceding are notintended to be exhaustive, but rather are among many which demonstratethat although significant attention has been devoted to increasing theinformation produced in analyzing time-series data, specificallyacoustic data, prior attempts do not satisfy the need for analysis ofdata across a wide range of frequencies.

BRIEF SUMMARY OF ONE EMBODIMENT

One embodiment which is intended to accomplish the foregoing objectscomprises a method for analyzing and interpreting time-series data. Themethod comprises acquiring data, computing the wavelet transform of thedata, computing the Lipschitz coefficients, that will provide thecharacterization of the singularities of the data, and representingthese singularities with respect to time. This method allows forinterpretation over recorded frequencies, also referred to asmulti-scale analysis. This type of analysis and interpretation isuseful, for example, in oil well borehole sonic data processing due tothe complexity of the data that is collected.

THE DRAWINGS

Objects and advantages of the present disclosure will become apparentfrom the following detailed description of embodiments taken inconjunction with the accompanying drawings, wherein:

FIG. 1 is a schematic view of the context in which one embodiment of thedisclosed system is intended to operate.

FIG. 2 is a flow chart of the steps of one embodiment.

FIG. 3 is a graphical representation of synthetic test data.

FIG. 4 is a graphical representation of local modulus maxima of datadepicted in FIG. 3.

FIG. 5 is a graphical representation of the Lipschitz Coefficients ofdata shown in FIG. 3.

FIG. 6 is a graphical representation of the accuracy of the estimationof the data shown in FIG. 3

FIG. 7 is a graphical representation of a set of synthetic sonic data.

FIG. 8 is a graphical representation of the singularities of the data inFIG. 7.

FIG. 9 is a graphical representation of the time shifted synthetic sonicdata of FIG. 7.

FIG. 10 is a graphical representation of the singularities of the timeshifted synthetic sonic data.

FIG. 11 is a graphical representation of actual acoustic data gatheredfrom a borehole.

FIG. 12 is a graphical representation of the singularities of the signalshown in FIG. 11.

FIG. 13 is a graphical representation of the singularities shown in FIG.12 after filtering.

FIG. 14 is a graphical representation of wavelet transformedcompressional signals and extrema lines in the case of a sand interval.

FIG. 15 is a graphical representation of extrema line fit, wavelet mapand singularity calculation in the case of compressional signalsrecorded in a sand interval.

FIG. 16 is a graphical representation of wavelet transformedcompressional signals and extrema lines in the case of a shale interval.

FIG. 17 is a graphical representation of extrema line fit, wavelet mapand singularity calculation in the case of compressional signalsrecorded in a shale interval.

FIG. 18 is a graphical representation of the comparison of monopolesingularities computed in the case of sand and shale intervals.

FIG. 19 is a schematic representation of the set up used for rockphysics experiments.

FIG. 20 is a graphical representation of the variation of singularityfor the shear wave as a function of stress.

FIG. 21 is a graphical representation of singularities comparisonbetween the two previous samples in the absence of glycol.

FIG. 22 is a graphical representation of singularities comparisonbetween both samples fully saturated by glycol.

FIG. 23 is a graphical representation of the comparison of computedsingularities, for both samples, saturated by brine and without glycol.

FIG. 24 is a graphical representation of the comparison of computedsingularities, for both samples, saturated by brine and glycol.

DETAILED DESCRIPTION

Referring now to the drawings and particularly to FIG. 1, there is showna schematic illustration of one operational context from whichsonic/acoustic data are collected in a borehole. FIG. 1 is a schematicrepresentation of an Array Sonic Tool (such as the Dipole Sonic ImagingTool (DSI) (Schlumberger™), Sonic Scanner (Schlumberger™), etc.)comprising a transmitter section 102 having a pair of upper and lowerdipole sources 103 arranged orthogonally in the radial plane and amonopole source 104. A sonic isolation joint 105 connects thetransmitter section 102 to a receiver section 106 that contains an arrayof one or more receivers. An electronics cartridge 107 is connected atthe top of the receiver section 106 and allows communication between thetool and a control unit 108 located at the surface via an electroniccable 109. With such a tool, it is possible to make both monopole anddipole measurements. This tool may be implemented for wireline loggingor for logging while drilling (LWD). In fact, time-series data collectedin any manner can be analyzed using this method. Other implementationssuch as Schlumberger's Sonic Scanner™ may also be used to collect thedata. Regardless of the collection platform, the analysis of the datawill proceed as described below.

Turning to FIG. 2, a flowchart is shown representing the high levelsteps in the instant method. This is an introduction to the method thatwill be further discussed in detail below. The first step, labeled 202,is acquiring time-series data to be analyzed. In one embodiment, thetime-series data are acoustic data collected using an instrument such asthe one shown in FIG. 1. The method is also useful on any type of timeseries data including, but not limited to, seismic data, resistivitydata, pressure data, and temperature data. The acoustic data arecollected from the borehole and may give insight to the surroundingfluid and solid formations. Once the data are acquired, wavelettransform of the data is computed in step 204. Wavelet transformessentially translates a signal in the time domain to the time-frequencyspace where the local regularity of the signal can be isolated andutilized. The result of wavelet transform is then used in thecomputation of the Lipschitz Coefficients in step 206. Thesecoefficients are indicative of the regularity of the signal and giveinsight into where the data changes in form. In fact it allows the studyand analysis of abrupt changes present in the analyzed signal or in itsderivatives. From the Lipschitz coefficients, the singularities of thesample data are determined in step 208. These singularities correspondto shapes of transitions included in the sample data and therefore arelinked to the medium crossed by the recorded signal. Knowing this, instep 210, the calculated singularities can be compared to knownsingularities that correspond to materials expected to be present in theborehole. Note that there are various ways to make this comparison. TheKolmogorov-Smirmov tests [W.H. Press et al, 1988] may be used, or theMahalanobis distance, for example, also may be computed. This known setof singularities is referred to as the control set and is determinedthrough experimentation such as, prior experience in analyzing boreholedata, theoretical or numerical modeling, or analysis of data fromlaboratory experiments. From this comparison, in step 212, properties ofthe material surrounding the tool 100 and the borehole are determined.In the instant case, these properties include the composition of boththe fluid and the solid formations surrounding the tool.

Step 202 is performed by tool 100, once this data are acquired, thesubsequent steps may take place downhole or at a processing center at ornear ground level, for example, at a processing center located at thesurface or at an underground location accessible from the surface. Thefirst computational step involves the computation of the wavelettransform of the acquired data. The wavelet transform translates asignal represented in the time domain into a time-frequency space wherethe local regularity can be characterized by decomposing the signal intoelementary wavelets. The major characteristic of the wavelets is thatall members of a given wavelet family Ψ_(s,x), are generated bytranslating and dilating a given initial wavelet ψ which is called “themother wavelet.” The wavelet transform of a function f(t) is defined by

$\begin{matrix}{{{Wf}\left( {s,x} \right)} = {\int_{- \infty}^{\infty}{{f(t)}\Psi_{s,x}^{*}{{t}.}}}} & (1)\end{matrix}$

where the Wf(s,x) is the coefficient of wavelet transform, s is thecoordinate on the scale axis, and x is the coordinate on the time axis.The notation ψ_(s,x)* represents the complex conjugate of the motherwavelet. In practice, the wavelet transform can also be expressed asconvolution between the function f(t) and the mother wavelet ψ(x/s).Additional information regarding wavelet transform may be found in, forexample, Grossmann and Morlet, 1984, and Mallat and Hwang, 1992, whichdiscuss the use of continuous wavelet transform to computesingularities.

The wavelet transform projects the signal on a family of translated anddilated basis functions and these functions can be represented in L¹,or, L² norm as

$\begin{matrix}{{{For}\mspace{14mu} L\; 1}{\psi_{s,x} = {\frac{1}{\sqrt{s}}{\psi \left( \frac{x - t}{s} \right)}}}} & (2) \\{{{For}\mspace{14mu} L\; 2}{\psi_{s,x} = {\frac{1}{s}{\psi \left( \frac{x - t}{s} \right)}}}} & (3)\end{matrix}$

The parameters x and s are the translation and scale parameters,respectively. The Fourier transform of the mother wavelet satisfies thefollowing conditions

$\begin{matrix}{{\int_{0}^{+ \infty}{\frac{{{\hat{\psi}(\omega)}}^{2}}{\omega}{\omega}}} = {{\int_{- \infty}^{0}{\frac{{{\hat{\psi}(\omega)}}^{2}}{\omega }{\omega}}} = {C_{\psi} < {+ \infty}}}} & (4)\end{matrix}$

The localization and characterization of singularities are best donewith continuous wavelet transform, which is translation invariant andallows focus on the sharp variations present in the signal or in thetime derivatives of the signal. This is an important property of wavelettransform allowing the characterization of the local regularity offunctions which is often measured with Lipschitz coefficients.

The singularities of a function are often characterized by its Lipschitzcoefficient values. As used herein, “singularities” refer to adiscontinuity in the analyzed signal or its derivatives. As shown inMallat and Hwang, the singularities can be expressed with the valuesthat are given by estimation of the Lipschitz coefficient. An example ofa function with singularity is the following canonical function

f(x)=β(x−x ₀)^(α) ,x>x ₀  (5)

where (x−x₀)^(α) can have a singularity at the time x₀. Note that thiscanonical function is used for better understanding. In fact most of thesingularities in a signal can be expressed in the form of the canonicalfunction at least near the singularities. This function is described bydifferent values of the parameters, x₀, β, α. The wavelet transform ofthis function can be expressed as follows

Wf(s,x)=βW[(x−x ₀)^(α)](s,x)  (6)

where (s,x) is the coordinate on the frequency and time axes. Here, theterm “frequency” has been used generally in place of a more precise“scale” in wavelet transforms.

By using the nth derivative of the mother wavelet ψ as the new motherwavelet and the L2 norm, the wavelet transform of Equation (6) can bedecomposed as follows:

$\begin{matrix}\begin{matrix}{{{\beta \left\lbrack {W\left( {x - x_{0}} \right)}^{\alpha} \right\rbrack}\left( {s,x} \right)} = {\beta \; {\frac{1}{s} \cdot \frac{^{n}}{\left( {x/s} \right)^{n}} \cdot {\psi \left( \frac{x}{s} \right)}}*\left( {x - x_{0}} \right)^{\alpha}}} \\{= {\beta \; {s^{\alpha} \cdot \frac{^{({n - \alpha - 1})}}{\left( {x/s} \right)^{({n - \alpha - 1})}} \cdot {\psi \left( \frac{x}{s} \right)}}*\frac{^{({\alpha + 1})}}{x^{({\alpha + 1})}}}} \\{\left( {x - x_{0}} \right)^{\alpha}} \\{= {{{\Gamma \left( {\alpha + 1} \right)} \cdot \beta}\; {s^{\alpha} \cdot {\psi^{n - \alpha - 1}\left( \frac{x}{s} \right)}}*{\delta \left( {x - x_{0}} \right)}}} \\{{= {{\Gamma \left( {\alpha + 1} \right)}\beta \; s^{\alpha}{\psi^{n - \alpha - 1}\left( \frac{x - x_{0}}{s} \right)}}},}\end{matrix} & (7)\end{matrix}$

where n≧α+1, the symbol * represents the convolution operator on thevariable x, ψ^(n−α−1) denotes the (n−a−1) derivative of ψ. In order toderive the previous expression, the following property has been used[Gel'fand and G. E. Shilov, 1962], i.e., for any function g(x) which iscontinuous and its 1st derivative is continuous,

$\begin{matrix}{{{g(x)}*\frac{^{\alpha + 1}}{x^{\alpha + 1}}\left( {x - x_{0}} \right)^{\alpha}} = {{\Gamma \left( {\alpha + 1} \right)}{\delta \left( {x - x_{0}} \right)}*{g(x)}}} & (8)\end{matrix}$

Once the precise location of a singularity at x₀ is known, the order ofthe Lipschitz regularity a can be estimated. The value α is also namedas Lipschitz coefficient. A function f(x) having a singularity at x₀with the order of Lipschitz regularity α is usually termed as “f(x) isLipschitz α at x₀.” The value of α for some simple functions withsingularities can be found based on Equations (7) and (8). For example,α=−1 for a dirac delta function δ(t), α=0 for a step function (orHeaviside function) H(t), α=1 for a ramp function r(t), α=2 for u(t),α=0.5 for a square root function sq(t), α=1.5 for sq3(t), and α=2.5 forsq5(t). Here, the definitions of these functions are:

${{\delta (t)} = {{0\mspace{14mu} {if}\mspace{14mu} t} \neq 0}},{{H(t)} = \left\{ {{{\begin{matrix}{{{+ 1}\mspace{14mu} {if}\mspace{14mu} t} > 0} \\{{\frac{1}{2}\mspace{14mu} {if}\mspace{14mu} t} = 0} \\{{{0\mspace{14mu} {if}\mspace{14mu} t} < 0},}\end{matrix}{r(t)}} = {\int{{H(t)}{t}}}},{{u(t)} = {\int{{r(t)}{t}}}},{{{sq}(t)} = t^{1/2}},{{{sq}_{3}(t)} = t^{3/2}},{{{sq}_{5}(t)} = {t^{5/2}.}}} \right.}$

The relationship between singularities and modulus maxima of the wavelettransform [Mallat and Hwang, 1992] is now discussed. Let n be a positiveinteger and α≦n. Let f(x)εL²(R). If f(x) is Lipschitz α at x₀, thenthere exists a constant A such that for all points x in a neighborhoodof x₀ and any scale s, such as

|Wf(s,x)|≦A(s ^(α) +|x−x ₀),  (9)

where Wf (s,x) represents the coefficients of wavelet transform.Conversely, let α<n be a non-integer value. The function f(x) isLipschitz α at x₀ if the following two conditions hold:

-   -   1) There exists some ε>0 and a constant A such that for all        points x in a neighborhood of x₀ and scale s

|Wf(s,x)|≦As ^(ε),  (10)

and

-   -   2) there exists a constant B such that for all points in a        neighborhood of x₀ and any scale s

$\begin{matrix}{{{{Wf}\left( {s,x} \right)}} \leq {{\beta \left( {s^{\alpha} + \frac{{x - x_{0}}}{{\log {{x - x_{0}}}}}} \right)}.}} & (11)\end{matrix}$

Knowing this, it is important to then characterize a particular class ofisolated singularities from the behavior of the wavelet transformmodulus maxima. Modulus maxima is the location of the local maxima ofthe absolute value of the wavelet transform Wf(s,x). Normally modulusmaxima coincide with the location of maximal energy. Let f(x) be afunction whose wavelet transform is defined over ]a,b[ and let x₀ε]a,b[.Assuming that a scale s₀>0 and a constant C exist such that for xε]a,b[and s<s₀, all the modulus maxima of Wf(s,x) belong to a cone defined by

|x−x ₀ |≧C·s.  (12)

Then, at all points x₁ε]a,b[,x₁≠x₀, f(x) is uniformly Lipschitz n in theneighborhood of x₁. Letting α<n be a non-integer, the function f(x) isLipschitz α at x₀ only if there exists a constant A such that at eachmodulus maxima (s,x) in the cone defined by

|Wf(s,x)≦As ^(α)  (13)

this theorem can also be written in an equivalent form

log|Wf(s,x)|≦log(A)+α·log(s).  (14)

The above equations are also valid for the case where α is anon-positive integer, for which the function also have singularity withLipschitz regularity α.

This form proves that the Lipschitz regularity (i.e., a) at x₀ is themaximum slope of a straight line that remains above log Wf|(s,x)|, on alogarithmic scale. This series of equations shows that the Lipschitzregularity is computed by finding the coefficient α that bestapproximates the decay of |Wf(s,x)| over a given range of scales largerthan 1. Once the mother wavelet is selected, the local modulus maximaare detected as it provides direct information on the regularity of thesignal as presented above. It is known that the modulus of the wavelettransform is biased in the presence of noise. Therefore, rather thanusing the standard approach consisting of the detection of the maxima ofWf(s,x), a more robust alternative approach is proposed and presentednow. In this case, the local regularities of the signal are estimatedthrough the computation of the sum of the modulus of its waveletcoefficients inside the corresponding “cone of influence” (COI), whichis the support of the wavelet at different scales, defined as follows:

$\begin{matrix}{{{S_{W}{f(x)}} = {\int_{{{x - x_{0}}} \leq {K \cdot s}}{{{{Wf}\left( {s,x} \right)}}{s}}}},} & (15)\end{matrix}$

where S_(w)f(x) is an operator equal to the sum of the modulus of thewavelet coefficients, |Wf(s,x)| considered in an interval |x−x₀|≦K·Salong the scale axis. K is the size of the interval which is taken asthe support of the mother wavelet considered in the analysis. Theposition and value of the local maxima of S_(w)f(x) will be used tocompute the Lipschitz computation rather than the local modulus maximaof the wavelet transform. The summation along the frequency axis isselected for the detection method as it provides a robust way to obtainthe estimation of the Lipschitz coefficient as the detection of themaxima will have a smaller perturbation caused by noise. By thistransformation, useful information in the signal is retained, and thenoise component in the signal is decreased. In consequence, the localmodulus maxima positions will be detected using this method havingsignificantly reduced noise perturbation. Using Equation 14 and the “newmaxima” approach described above, estimation of the Lipschitzcoefficients may be obtained in a robust manner.

An illustrative example of the method can be shown by considering asynthetic signal built as a linear combination of signals defined usingthe canonical expression presented in Equation 5. The signal is composedof components exhibiting various singularities order at various timepositions. The signal test is defined by the following equation:

s(t)=δ(t−t ₁)+0.5H(t−t ₂)+0.01r(t−t ₃)−0.001sq ₃(t−t ₄)+0.0001u(t−t₅)−0.00001sq ₅(t−t ₆)  (16)

where δ(t), H(t), r(t), u(t), sq(t), sq₃(t), and sq₅(t). have beendefined previously. The test signal contains various singularities

$\alpha = \left\{ {{- 1},0,1,\frac{3}{2},2,\frac{5}{2}} \right\}$

located at the time t={t₁, t₂, t₃, t₄, t₅, t₆}={255, 511, 767, 1023,1279, 1379}. A graphical representation of the signal is shown in FIG.3.

FIG. 4 shows the modulus computed for the signal. Here the horizontalaxis is time shift, x, and the vertical axis is scale, s. The lines ofmodulus maxima are then obtained and used to compute the order of thesingularity. Note that the maxima lines will converge to the location ofthe singularity (e.g., x0 in Equations 5 and 6) as the scale, s, issmall. Accordingly, this approach could be useful to detect theexistence of various components of a signal and their respectiveproperties related to the material being investigated.

A log-log plot representation gives a direct estimation of Lipschitzcoefficient α using linear regression. The coefficient α will correspondto the slope of this linear fit. FIG. 5 shows the estimation of theLipschitz coefficient, α, on the vertical axis and its localizationtime, x, on the horizontal axis. FIG. 6 shows the accuracy of theestimation for the synthetic signal. Note that all the varioussingularities have been properly located in the time domain and theirrespective order properly estimated.

Singularities have been used in the past for various other applications,like image processing, for example, as discussed in X. H. Wang et al, K.A. Innanen, and C. L. Du and W. L. Hwang. Disclosed herein is a newapproach to analyze and extract singularities as attributes from sonicdata and other time-series data. The Lipschitz coefficients α obtainedcorrespond to the shapes of transitions in the original waveform andtherefore are indicative of the medium encountered by the acoustic wavesdetected by the receiver(s) and recorded by the recording system. Thelocation of a in time indicates that at that time, the waves propagatingin the formation or fluid changed or encountered a particular medium.Conventionally time series data are represented by the amplitudevariation as a function of time. The extraction of Lipschitzcoefficients α and the singularity locations is a way to quantifyvariations due to the propagation in the medium. As it is possible toquantify and detect various components, it is therefore possible todenoise and filter the time series of interest using the detectedsingularities.

One method for filtering is now discussed. One comparatively simpleapproach uses a predefined threshold. Singularities that have a valuegreater than or less than the predetermined threshold value areretained. The predefined threshold may be set by the user, or it may bedefined based on experience, i.e., after analyzing a certain number ofdata sets it is possible to define a certain value for certain types ofrocks. Another approach is based on automatic filtering, without humanintervention. In this approach, all the computed singularities in onesection of a well are used and a clustering analysis is applied to thedata. The clustering analysis will extract various classes ofsingularities that are present in the analyzed section, i.e., thevarious components present in the data. Therefore, the algorithm willfilter the data corresponding to the various classes found by theclustering analysis.

It can be shown that a shift in time or a change in amplitude will notchange the results of Lipschitz coefficients α. FIG. 7 is arepresentation of a synthetic sonic signal. Using the previouslydescribed method, the singularities shown in FIG. 8 are detected. FIG. 9shows the signal after it has been time shifted and FIG. 10 shows thesingularities determined after this shifting. A similar result isobtained when the signal is changed only in amplitude. Thissingularities detection demonstrates that the Lipschitz coefficients αcapture only the shapes of the waveform, neither shifting on the timeaxis nor the variation of amplitude have influence on α.

The application of the method to real borehole data is discussedhereinafter. FIG. 11 shows a set of acoustic data gathered from aborehole. By analyzing the singularities in the waveforms, it ispossible to filter out unwanted wave components. For example, the datashown in FIG. 11 represent at least two signals, the desired Ydipole,contaminated by a monopole signal. The singularities of the data aredetermined using the method described above. These singularities arerepresented in FIG. 12 as a function of time and depth. As can be seen,the singularities track the shape of the original signal. Further, bylooking at the variety of values, it becomes clear that what appeared tobe each waveform is actually a combination of waveforms that may be theresults of interference by undesirable waves. One embodiment of thisinvention is to remove the undesirable interferences by removing theextraneous singularities from the original waveforms. FIG. 13 shows thesingularities after the removal of spurious singularities 501 and 502.This filtering process can be very beneficial to improve the quality ofthe measured data.

A second example is presented to illustrate another application of theproposed method. Two sets of monopole waveforms recorded with anacoustic tool as described above have been acquired respectively in ashale and a sand section of a well. First, the case of monopole recordedin sand is presented. FIG. 14 presents the compressional signal(middle), its wavelet transform (top), and the extracted maxima lines(bottom) in the case of the sand interval. FIG. 15 presents how thesingularity is computed, i.e., the fit in the log|Wf(s,x)| versus log(scales) domain. FIGS. 16 and 17 present similar results, but in thecase where the compressional wave has been recorded in the shalesection. Note that there are visible differences between both figuresillustrating the singularity difference values between both waveforms.FIG. 18 presents the singularity comparisons for compressional wavesrecorded in a shale and sand interval. Note the clear difference betweenthe singularities computed for both cases, i.e., it is possible todifferentiate a waveform recorded in a sand section of the well from onerecorded in a shale section based on its singularity values. Bothexamples illustrate clearly the use of singularity to discriminatequantitatively waveforms recorded in different environments, i.e.,recorded in a substratum with different petrophysical properties.

In the following examples, singularities in waveforms are shown todiscriminate various rock properties, such as: stress effects, rocktypes, pore fluid types, etc. In this case, the results for S waves aredemonstrated to illustrate the ability to evaluate rock properties byuse of the instant method. Similar results can also be obtained with Pwaves. Ultrasonic measurements were performed with an acousticpressurized cell designed by Temco Inc [Temco]. FIG. 19 shows the setupwith the position of transducers [NER] that enable theemission/reception of ultrasonic waves (P-, S₁- and S₂-waves that areorthogonally polarized) in a frequency range of 1 kHz to 4 MHz. Thepressure control is made with three separated pumps [ISCO] for the porepressure P_(p), the radial (or confining) pressure P_(r), and the axialpressure P_(a), respectively. The sample is isolated inside a Hasslersleeve to avoid any contact with the external fluid allowing thepressure control. As one example, 2 different carbonates rocks wereselected because of their wide range of permeability (Table 1 below).

Porosity Grain density Permeability Sample Φ (%) ρ_(g) (g/cc) k (mD)Q244 19.77 2.720 1.26 Q263 18.77 2.901 343.72

The Q244 presents quasi-spherical pores and is supposedly non-connectedregarding its permeability. The Q263 is composed of a wide range of poreshapes, which are connected together.

The first analysis concerns the study of the sample Q244 as a functionof stress without fluid inclusion in the system (0% Glycol). In thiscase, only the first singularity corresponding to the S wave is plotted.It can be clearly observed that FIG. 20 shows the variation ofsingularity with the applied stress. Note the singularity valuesincrease to reach a plateau at higher stress.

Next the two samples presented previously are compared using theirsingularity values. In this case, the shear wave recorded for eachsample is used. In absence of glycol, i.e., no fluid has been injectedin the rock, FIG. 21 presents the computed singularities for bothsamples under an applied stress of 2000 and 3000 psi, respectively. Notethat the singularities of the Q263 sample have higher values than theQ244 sample. Next, the case where the rock is 100% saturated by glycolis considered. The same applied stress is set. In both cases thesingularities of the Q263 sample have higher values than the Q244sample, but they are smaller than in the previous case with no fluid(note FIG. 22).

A final test consists in saturating the rock with brine and comparingthe singularity calculations with the ones obtained when the rock is dryor saturated with glycol. In this example, the influence of fluid typeon the singularities values is studied. In this case, two pressures areconsidered (2000 and 3000 psi). FIGS. 23 and 24 clearly show thevariation of singularities depending on the presence or absence of fluidand on the fluid content. It has been shown in this case that it ispossible to use singularities to discriminate rocks properties as afunction of fluid type, rock type, permeability, porosity, stress, etc.

For further analysis, it may be desired to filter out the noise in thesignal; the use of singularities makes this possible because certainsingularities known to be associated with noise can be removed. Themajority of the singularities in regions 501 and 502 (note FIG. 12) arefiltered and removed from the analysis. The results of this removal areshown in FIG. 13. This keeps only the singularities that are of interestin the determination of the materials surrounding the receiver. Thesingularity analysis of an acoustic signal is a powerful tool inanalyzing the properties of formations surrounding the borehole.Knowledge of the formation properties, i.e., for answer products, isessential to understand and exploit the reservoir. An answer product maycontain a wide range of information including, but not limited to,reservoir parameters, seismic attributes, reservoir mapping, scaling,formation type, level of heterogeneity of surrounding rock, porosity,rock fraction detection, FMI imaging, and image processing. Furthermore,by adding other known petrophysical measurement data, such as data thatare typically acquired through various tool measurements, for example,resistivity, nuclear, density, permeability, electro-magnetic, and thelike, to the singularity analysis additional information about theformations may be provided. For example, by adding data from resistivitymeasurements with the data from singularities will increase the level ofinterpretation of the reservoir. As previously mentioned, various typeof data known to those skilled in the art may be utilized for thepurposes described herein.

An additional use for the algorithm is in extracting information ofinterest from the determined singularities. This is similar to the abovewhere an individual signal is separated from a group of other signals ofinterest. Additionally, patterns of interest in images can be discernedwhich is useful in a wide range of applications. Once a signal or imageis collected, it is also possible to characterize a discontinuity ofinterest and determine when the formation or other material changes.

In the foregoing description, reference has been made to a variety ofembodiments. The subject disclosure, however, is not limited to wellholetechnology and is rather intended to provide useful application in allcontexts where time series data analysis is desired. Those skilled inthe art and familiar with the instant disclosure may recognizeadditions, deletions, modifications, substitutions, and other changeswhich fall within the purview of the disclosure and claims.

1. A method for analysis of oil field time-series data, comprising:acquiring data in an oil field; extracting singularities of the data;using singularities to interpret the formation properties related to thedata.
 2. The method as defined in claim 1, wherein the acquired datacomprises one or more of sonic data, seismic data, and ultrasonic data.3. The method as defined in claim 2, wherein the seismic data comprisesvertical seismic profile (VSP) data.
 4. The method as defined in claim1, wherein the formation properties comprise one or more of rock type,fluid type, permeability, and porosity.
 5. The method as defined inclaim 1, further comprising the steps of: comparing said singularitiesto a control set; and determining said formation properties based onsaid comparison.
 6. The method as defined in claim 1, further comprisingthe step of: filtering based on said singularities.
 7. The method asdefined in claim 1, wherein: prior to said extracting step: computing awavelet transform with said acquired data; and computing Lipschitzcoefficients utilizing said wavelet transform.
 8. The method as definedin claim 7, wherein: computing said Lipschitz coefficients comprisesdetermining local modulus maxima; and said local modulus maxima arecomputed by summing along the frequency axis.
 9. The method as definedin claim 8, wherein: said summing of local modulus maxima are computedusing the following formula:S_(W)f(x) = ∫_(x − x₀ ≤ K ⋅ s)Wf(s, x)s.
 10. The method asdefined in claim 1, wherein: said data is acquired utilizing at leastone source and at least one receiver.
 11. The method as defined in claim1, wherein said time series data is collected at or above ground level.12. The method as defined in claim 1, wherein: said time series data iscollected downhole in a borehole.
 13. The method as defined in claim 1,wherein: said time-series data is one or more data relating toresistivity, nuclear, density, permeability, electro-magneticmeasurements.
 14. The method as defined in claim 1, wherein: saidtime-series data is acoustic data.
 15. The method as defined in claim 1,wherein: said time-series data is petrophysical data acquired throughoil field tool measurements.
 16. The method as defined in claim 15,wherein said time series data is collected via wireline.
 17. The methodas defined in claim 15, wherein said time series data is collected whiledrilling.
 18. The method as defined in claim 1, wherein saidinterpretation is performed downhole.
 19. The method as defined in claim1, wherein said interpretation is performed at or near ground level. 20.A method for analysis of oil field time-series data comprising:acquiring data in an oil field; extracting singularities of the data;identifying the singularities related to unwanted signals; filtering outunwanted signal related to the identified singularities.
 21. A methodfor analysis of oil field time-series data comprising: acquiring data inan oil field; extracting singularities related to the signals; usingsingularities in combination with other measurements to interpret theformation properties.